https://bioinformatics-core-shared-training.github.io/cruk-summer-school-2018/RNASeq2018/html/06_Gene_set_testing.nb.html

# Load library and search for the organism
library(clusterProfiler)
library(pathview)
library(tidyverse)
search_kegg_organism('sasa', by='kegg_code') ## Salmo salar = sasa in the KEGG data base

PATHWAY ANALYSIS FOR POMV 6 HPI

The following code loads ALL differentially expressed genes with an adjusted p-value < 0.05, and then selects those with an absolute logfold change of 2.

POMV6 <- read.csv('Results/ControlvsPOMV6_ALL.csv') 
POMV6_SigGeneList <- POMV6[abs(POMV6$logFC) >= 2,]
POMV6_SigGeneList <- POMV6_SigGeneList[order(abs(POMV6_SigGeneList$logFC), decreasing = TRUE),]
rownames(POMV6_SigGeneList) <- NULL
POMV6_SigGeneList_logFC <- POMV6_SigGeneList$logFC
names(POMV6_SigGeneList_logFC) <- POMV6_SigGeneList$ENTREZID
kk_POMV6 <- enrichKEGG(gene = names(POMV6_SigGeneList_logFC), organism = 'sasa')
head(kk_POMV6)
barplot(kk_POMV6)

dotplot(kk_POMV6)

emapplot(kk_POMV6)

cnetplot(kk_POMV6, categorySize="genNUM", foldChange=POMV6_SigGeneList_logFC, circular = TRUE)

PATHWAY ANALYSIS FOR POMV 24 HPI

POMV24 <- read.csv('Results/ControlvsPOMV24_ALL.csv') # CHECK FOR NAs
POMV24_SigGeneList <- POMV24[abs(POMV24$logFC) >= 2,]
POMV24_SigGeneList <- POMV24_SigGeneList[order(abs(POMV24_SigGeneList$logFC), decreasing = TRUE),]
rownames(POMV24_SigGeneList) <- NULL
POMV24_SigGeneList_logFC <- POMV24_SigGeneList$logFC
names(POMV24_SigGeneList_logFC) <- POMV24_SigGeneList$ENTREZID

kk_POMV24 <- enrichKEGG(gene = names(POMV24_SigGeneList_logFC), organism = 'sasa')
head(kk_POMV24)
barplot(kk_POMV24, showCategory=20)

dotplot(kk_POMV24, showCategory=20)

emapplot(kk_POMV24)

In the following plot I removed node labels because they were too many …

library(igraph)
g <- cnetplot(kk_POMV24, categorySize="genNUM", foldChange=POMV24_SigGeneList_logFC, layout = 'nicely', node_label = FALSE)
plot(g, vertex.label.cex=0.5)

PATHWAY ANALYSIS FOR ISAV 6 HPI

ISAV6 <- read.csv('Results/ControlvsISAV6_ALL.csv') # CHECK FOR NAs
ISAV6_SigGeneList <- ISAV6[abs(ISAV6$logFC) >= 1.5,]
ISAV6_SigGeneList <- ISAV6_SigGeneList[order(abs(ISAV6_SigGeneList$logFC), decreasing = TRUE),]
rownames(ISAV6_SigGeneList) <- NULL
ISAV6_SigGeneList_logFC <- ISAV6_SigGeneList$logFC
names(ISAV6_SigGeneList_logFC) <- ISAV6_SigGeneList$ENTREZID

kk_ISAV6 <- enrichKEGG(gene = names(ISAV6_SigGeneList_logFC), organism = 'sasa')
head(kk_ISAV6)
barplot(kk_ISAV6)

dotplot(kk_ISAV6)

emapplot(kk_ISAV6)

cnetplot(kk_ISAV6, categorySize="genNUM", foldChange=ISAV6_SigGeneList_logFC)

PATHWAY ANALYSIS FOR POMV 24 HPI

# Loaded gene lists como from Dif_expression (ALL differentially expressed genes with an adjusted p-value < 0.05)
ISAV24 <- read.csv('Results/ControlvsISAV24_ALL.csv') # CHECK FOR NAs
ISAV24_SigGeneList <- ISAV24[abs(ISAV24$logFC) >= 1.5,]
ISAV24_SigGeneList <- ISAV24_SigGeneList[order(abs(ISAV24_SigGeneList$logFC), decreasing = TRUE),]
rownames(ISAV24_SigGeneList) <- NULL
ISAV24_SigGeneList_logFC <- ISAV24_SigGeneList$logFC
names(ISAV24_SigGeneList_logFC) <- ISAV24_SigGeneList$ENTREZID

kk_ISAV24 <- enrichKEGG(gene = names(ISAV24_SigGeneList_logFC), organism = 'sasa')
head(kk_ISAV24)
barplot(kk_ISAV24, showCategory=20) ## SHOW ALL CATEGORIES!

dotplot(kk_ISAV24, showCategory=20)

emapplot(kk_ISAV24, showCategory=20)

cnetplot(kk_ISAV24, categorySize="genNUM", foldChange=ISAV24_SigGeneList_logFC)